Particle-in-Cell simulations of electron spin effects in plasmas 
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We have here developed a particle-in-cell code accounting for the magnetic dipole force and for the magne- 
tization currents associated with the electron spin. The electrons is divided into spin-up and spin-down pop- 
ulations relative to the magnetic field, where the magnetic dipole force acts in opposite directions for the two 
species. To validate the code, we have studied the wakefield generation by an electromagnetic pulse propagating 
parallel to an external magnetic field. The properties of the generated wakefield is shown to be in good quan- 
titative agreement with previous theoretical results. Generalizations of the code to account for more quantum 
effects is discussed 
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Recently much work has been devoted to the field of quan- 
tum plasmas. Many aspects of the field are well described 
in books and review papers [1-5]. The interest in the field 
has been motivated by laboratory applications in for example 
plasmonics [6, 7], quantum wells [8] and high -density laser 
plasma interaction [3, 9] as encountered e.g. in inertial con- 
finement fusion (ICF) schemes. Furthermore, astrophysical 
applications [10] of quantum plasmas has played a role. The 
theoretical approaches includes quantum kinetic theories such 
as e.g. the Kadanoff-Baym kinetic [1] or the Wigner equation 
[2], various quantum hydrodynamic equations [11-13] and 
spin kinetic theories [14-17]. Due to the complex nature of 
the governing equations, many problems involving nonlinear- 
ities and/or inhomogeneities must be addressed numerically. 
A successful method applied to classical plasmas is the cele- 
brated particle-in-cell (PIC) approach (see e.g. [18-21]). Due 
to the classicality of the concepts this have had a limited im- 
pact to quantum plasmas. It should be noted, however, that the 
Feynman path integral formulation has been used to develop a 
PIC treatment [22, 23] that includes particle dispersive effects 
in a semi-classical fashion. This is computationally costly, 
and the number of quantum particles that can be included in 
the code is far less than in the classical case [22, 23]. 

In our paper we have chosen another approach in order to 
develop a PIC-code including certain quantum features. We 
first note that the Wigner equation reduces to the classical 
Vlasov equation for macroscopic scale lengths longer than 
the thermal de Broglie wavelength. Hence particle disper- 
sive effects [2, 24] can be neglected on such scales. On the 
other hand, the physics associated with the electron spin (e.g. 
the magnetic dipole force and the spin magnetization currents 
[15-17]) does not vanish for long scale lengths. Hence, in 
this paper we will aim to develop a PIC code applicable on 
macroscopic scales, accounting for the magnetic dipole force 
and the spin magnetization currents. While a classical mag- 
netic dipole moment fits very well into the PIC-concept, it is 
clear that the difference between a classical magnetic dipole 
moment and the spin must be acknowledged. The present ap- 
proach builds on the findings of Ref. [25]. It was then ob- 
served that more elaborate models for the spin physics reduces 
to a simple one for frequencies below the spin precession fre- 
quency (which is approximately the same as the cyclotron fre- 



quency). Specifically this meant that for a dynamical time 
scale slower than the cyclotron frequency, the electrons could 
be modelled as consisting of two fluids, one with a spin up 
state relative the magnetic field and one with a spin down 
state. The two electron species are then subject to magnetic 
dipole forces acting in opposite directions. A ID PIC-code 
based on this concept has been developed, and is tested in 
the present paper. In particular we let a short electromagnetic 
(EM) pulse propagate along an external magnetic field in a 
plasma, and study the wake field induced. Comparing with 
theoretical results, we find that the PIC-code is able to repro- 
duce previous findings [26, 27]. 

Due to the classical nature of the PIC-concept, it is theo- 
retically difficult to generalize it to the quantum regime. Al- 
though the basic effect of wave particle dispersion (as de- 
scribed by the Schroedinger equation) can be handled to a 
certain extent (see e.g. Refs. [22, 23]) it is computationally 
costly. In particular the number of particles needed to capture 
particle dispersion is increased by a factor of at least a thou- 
sand. On the other hand, the effect of particle dispersion is not 
important for scale lengths much longer than the characteris- 
tic de Broglie wavelength [2, 16]. However, other quantum 
effects may play a role for the long scale dynamics [15-17]. 
In particular, when HbB/IcbT or /igB/zn/c^ approaches unity 
(here /ig = eh/2m e is the Bohr magneton, B is a characteris- 
tic magnetic field value, T is a characteristic temperature, e 
is the elementary charge, h is the Planck constant divided by 
2n, m e and m,- is the electron and ion mass respectively, and 
ca is the characteristic Alfven velocity), the magnetic dipole 
force and the magnetization current associated with the elec- 
tron spin may become important [15, 28]. To a significant 
extent, the magnetic dipole force associated with the spin is 
comparable to a classical magnetic dipole force, and the mag- 
netic moment associated with the magnetic field is compara- 
ble to a classical magnetic moment. This means that spin ef- 
fects can be conceptually straightforward to include in a PIC- 
scheme. A complication is that unlike a classical magnetic 
dipole moment, the probability distribution of the spin for a 
single particle is by necessity spread out (since the different 
spin components does not commute). This suggests that spin- 
velocity correlations [29] may affect the dynamics, and that 
updating the variables using the Heisenberg equation of mo- 
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tion may be somewhat inaccurate. However, the case where 
the dynamical time scale is slower than the electron spin pre- 
cession frequency (which is approximately the same as the 
electron cyclotron frequency) and the spatial scale is longer 
than the Larmor radius has been studied by Ref. [25]. It was 
then deduced that the electrons can be described by a rather 
simple model, where only two spin states, up and down rel- 
ative to the magnetic field are needed. The quantum effects 
then enter as a magnetic dipole force in the equation of mo- 
tion (in opposite directions for the two spin states, which is 
described as two different species), and as a magnetization 
along the magnetic field, which is proportional to the density 
difference between the spin-up and spin-down electrons. 

Based on the concepts described above, we will extend an 
ID Particle-In-Cell simulation (LPIC++) (see Ref. [21]) in 
order to incorporate the effect of up- and down spin states in 
the governing dynamics. In this model the ions are treated 
as a classical particles. However the electrons come in two 
different species, one with spin up and the other with spin 
down. The presence of an external magnetic field along the 
direction of laser propagation is also included in the model. 
The governing equations for particle motion with two electron 
species will now read as, 



d\ i 

m e ~r = 9«(E + x B) - jU B VB 
dt 



d\f 

m e —^- = q e (E + \^ x B)+/1 B VB 
dt 



(1) 



(2) 



where the arrows indicate the different spin states and B is the 
magnitude of the magnetic field. Moreover Ampere's law is 
written as, 



1 <9E 

Vx(B- J uoM) = J uoj/ + ^-^- 



(3) 



where the free current including the ions is jf = q e [n±Y± + 
«|V|) + qj\i, and the magnetization due to spin is M = 
/Xflb(nj — rtf), where b is the unit vector in the direction of 
magnetic field, i.e. b = B/B. Two assumptions is necessary 
for this model to hold 

i) That the spin-relaxation time is longer than the time 
scales of study. Otherwise the two spin-states cannot be 
treated as independent fluids. 

ii) That the time scale of study is longer than the spin pre- 
cession time scale. Otherwise the spin vector (and the mag- 
netic moment) will not necessarily point in the direction of the 
magnetic field. 

Before we start the studies it is convenient to rewrite the 
equations in dimensionless form. We use the inverse laser fre- 
quency o -1 to normalize time (i.e. the normalized time t„ 
is t n = cot), c/ co to normalize distance, the normalized spin- 
down (up) density is n n u^) = (n^\e 2 / £om) j CO 1 , and the di- 
mensionless electric and magnetic field are E„ = eE/m e COc 
and B„ = eB/m e (0, respectively. Dividing the magnetic field 
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FIG. 1. Simulation geometry of the problem. 



into a homogeneous static field Bo = Bqx and a dynamic mag- 
netic field, the basic equations in normalized form becomes 



dt„ 



/ co c ^\ hco „ 
B -x 7 VB 

\ co J mcr 



dt n 
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and 



V x ( B - -^b(n; -«t) ] = (rtjv; - 



«t v t) + -57 ( 6 ) 



Here for notational convenience the subscript n denoting nor- 
malized quantities have been omitted on all variables and op- 
erators. The term (co c / ca)x, where co c = eBo/m represents the 
constant static part of the magnetic field, and hence B is the 
dynamic part only. Consequently we have B = |(fi) c /fi))x + B|. 
For all simulations presented below the dimensionless param- 
eter hco/mc 2 is chosen to be hco/mc 2 = 2k x 10~ 4 . The typ- 
ical simulation geometry for the problem is shown in Fig. 
1. For all the results presented here the normalized electric 
field amplitude is considered as 0.15. The circularly polar- 
ized laser with FWHM duration of 10 period times propagates 
through the plasma along the x direction. The external mag- 
netic field is considered along the x direction with amplitude 
in dimensionless unit (bo = eBo / m e co). With these parameters 
the characteristic velocities of the problem (e.g. phase veloc- 
ity and group velocity) will be of the order of c.The thermal 
velocity is considered to be much smaller than this, and hence 
the temperature is effectively zero in our simulations. Finally, 
the plasma of length 80 wavelengths is taken to have the total 
unperturbed electron density of + n-j- = 0.3 (Fig. 1). 

In Fig. 2 the density perturbation of the electrons with 
spin-up and with spin-down is shown after the passage of 
a circularly polarized EM wave. As can be seen the den- 
sity perturbations of the two species of electrons are almost 
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FIG. 2. Density perturbation with respect to unperturbed density for 
electron species with spin + 1/2 and —1/2 are presented at 60t with 
external magnetic field bo = 70. Inset shows the magnified version 
of the +1/2 spin state plot, where the presence of the small second 
harmonic component can be seen. 
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FIG. 4. Electron density response to the electric field of circularly 
polarized laser pulse is presented. Both electric field and difference 
between electron densities of both species are normalized to unity. It 
should be noted that both curve are identical and it is not possible to 
distinguish one from another in this figure. 



opposite. This is also clear from Fig. 3, where the total 
(spin-up + spin-down) density perturbation is shown, which 
is smaller than the individual density perturbations of the up- 
and down- species by roughly a factor 20. Furthermore, we 
see that the total density perturbation has a rapidly oscillating 
second harmonic component which is of the same order of 
magnitude as the low-frequency component. Since the elec- 
tromagnetic wave is circularly polarized, we note that the sec- 
ond harmonnic components are associated with the discrete 
particle (collisional) effects of the PIC-code. To understand 
the results in more detail, we first note that the spin proper- 
ties in the given geometry only enters nonlinearly. Since the 
magnetic field perturbation is purely transverse, we note that 

VB = vJbY+b\ « (1/2)V(B^/B )» including only up to 
second order nonlinearities in the transverse wave magnetic 
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FIG. 3. Spatial profile of total electron density at 60t for the external 
magnetic field &o = 70. 



field B±. The spin dependent part of the perturbation is then 
induced by the magnetic dipole force, which is directed in op- 
posite directions for the spin-up and down species. 

In order to verify the interpretation that the density differ- 
ence between the up- and down electron species comes is a 
direct response to the magnetic dipole force of the wave field, 
we compare the profile of |E (which is proportional to B ± ) 
with that of the density difference n+ — n~ . The result is 
shown in Fig. 4. It turns out that the matching is so accu- 
rate that two separate curves cannot be recognized, as they are 
fully overlapping. 

Next our aim is to make a quantitative comparison of the 
output from the PIC-code with analytical results derived by 
Ref. [27]. Eqs. (19) and (20) together with Eq. (20)ofRef. 
[27] show that for strong magnetic fields (co -C CO c ) the peak 
electric amplitude of the wakefield scales as 

ci(« 0+ -«o-) c 2 
E ^ Bo + ¥ Q 

where c\ and c 2 are constants. Fitting the peak amplitude as 
a function of the magnetic field as calculated from the PIC- 
code, with a curve of this type, a good agreement is found, 
as shown in Fig. 5. Specifically three different cases are in- 
vestigated. A pure spin up state, a pure spin down state and 
a mixed state with «o+ = «o--The constant c\ varies less than 
0.2% between the three cases, whereas c 2 varies slightly from 
5.735 x 10~ 4 to 5.328 x 10 -4 when switching from a pure 
spin-up plasma to a spin down plasma. 

Another interesting variable to study is the spin polarization 
induced by the EM wave for an plasma with equal number 
of particles in the spin up and the spin down states initially. 
The spin polarization is monitored by plotting the peak of the 
density difference n + — which is also closely related to the 
spin induced magnetization. According to the theory (see e.g. 
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FIG. 5. Maximum amplitude of wakefield at 60t is presented by 
varying the magnitude of external magnetic for electrons species 
with different spin states. For all the plasma electrons with spin 
down states (red circles) the associated function /(£>o) is fitted (blue 
line) with a = 5.328 x 1(T 4 and S, = 7.629 x 1(T 6 . For all the 
plasma electrons with spin up states (green circles) the associated 
function g(bo) is fitted (magenta line) with rj = 5.753 x 10~ 4 and 
<5 = — 7 . 6 1 5 x 1 CP 6 . For the mixed electron species present in plasma 
(black line) nicely coincides with the classical results (black open cir- 
cles) with fls = 0. 

Ref. [27]) the scaling with the magnetic field for (0 <C (O c is 

n+ — ti- oc Bq 
which is confirmed in Fig. 6 to a very good accuracy. 
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FIG. 6. The magnetization as a function of applied magnetic field 
(blue circles) is fitted with function /(2>o) (red line). 

Extensions of PIC-codes to include the spin effects can be 
particularly relevant for ICF plasmas. As demonstrated by 
Ref. [26], the spin dependent part of the ponderomotive force 
can give raise to significant spin polarization for plasmas with 
metallic densities and higher, provided the pulse length is not 
too short (cf Eq. (24) of Ref. [26]) Furthermore, the genera- 
tion of strong quasi-static magnetic field in laser-plasma inter- 
action, even approaching giga-gauss values [30], means that 



the background plasma can be significantly spin-polarized, 
since [IbB/IcsT can approach unity in the fast ignition sce- 
nario where the compression takes place before the heating 
[31]. In the present study we have focused on validating val- 
idate the code against the theoretical results of Ref. [27], and 
hence picked parameter values corresponding to astrophysical 
regimes. By comparing with analytical theory we have veri- 
fied that the code executes properly and with good accuracy. 
In particular, studying the effects of a weakly modulated high- 
frequency pulse, it has been shown that the low-frequency 
wake-field is modified by the magnetic dipole force as pre- 
dicted by theory (in particular, see Fig. 5). Furthermore, the 
spin polarization induced by the high-frequency wave in an 
system that initially is not spin-polarized agrees well with the- 
ory (see Fig. 6). This paper should be considered as a first step 
towards a full spin PIC-code. Besides generalizing the code to 
higher dimensions, it would be of much interest to include the 
spin dynamics (e.g. spin precession), in which case it would 
be possible to resolve the dynamics on the high-frequency cy- 
clotron scale. 
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